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Abstract. This work is concerned with the ill-posed inverse problem of estimating turbulent 
flows from the observation of an image sequence. From a Bayesian perspective, a divergence-free 
isotropic fractional Brownian motion (fBm) is chosen as a prior model for instantaneous turbulent 
velocity fields. This self-similar prior characterizes accurately second-order statistics of velocity 
fields in incompressible isotropic turbulence. Nevertheless, the associated maximum a posteriori 
involves a fractional Laplacian operator which is delicate to implement in practice. To deal with 
this issue, we propose to decompose the divergent-free fBm on well-chosen wavelet bases. As a 
first alternative, we propose to design wavelets as whitening filters. We show that these filters are 
fractional Laplacian wavelets composed with the Leray projector. As a second alternative, we use 
a divergence-free wavelet basis, which takes implicitly into account the incompressibility constraint 
arising from physics. Although the latter decomposition involves correlated wavelet coefficients, we 
are able to handle this dependence in practice. Based on these two wavelet decompositions, we 
finally provide effective and efficient algorithms to approach the maximum a posteriori. An intensive 
numerical evaluation proves the relevance of the proposed wavelet-based self-similar priors. 
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1. Introduction. This work is concerned with the ill-posed inverse problem of 
estimating turbulent motions from the observation of an image sequence. Turbulence 
motion phenomena are often studied in the context of incompressible fluids, which is 
the setting of this paper. This inverse problem arises in the context of experimental 
physical settings, where one is interested in recovering the kinematical state of an 
incompressible turbulent fluid flow from the observation of a sequence of images, e.g. 
particle image velocimetry in experimental fluid mechanics, wind or ocean currents 
retrieval from satellite imagery in geophysics. Solving accurately this type of inverse 
problems constitute an important issue since a complete physical theory is still missing 
for turbulence phenomenology. 

More specifically, the above inverse problem can be viewed as a Maximum A 
Posteriori (MAP) estimation of a vectorial field u over a space of admissible solutions: 

u*(yo,yi) e argmaxp5j,|uPu, 

where j/o ^md yi are two consecutive observed images, u is a velocity field, and 5y 
is a function of y^, yi and u, which characterizes the evolution between yo and yi. 
In this problem, u is not observed and the incompressibility constraint demands the 
motion field u to be divergent-free. In this Bayesian framework, psy\u denotes the 
likelihood model, which relates the motion of the physical system to the spatial and 
temporal variations of the image intensity. The adjunction of a prior information pu 
for the velocity field u is in this case mandatory since, as we will see in section [4] this 
non-linear problem is under-constrained. 
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Concerning the choice of psy\u: relevant physical models have been proposed for 
fluid flow imagery. A review in the context of experimental fluid mechanics can be 
found in [30|. We will assume in this paper the simple model where yo and t/i are 
the solutions between two consecutive times of a transport equation driven by u, see 
section |4] for more details. 

Concerning the choice of pu, the incompressible Navier-Stokes equations perfectly 
describe the structure of an incompressible velocity field, i.e. the prior for u. However, 
this implicit choice of constraints leads to an optimization problem which is often 
severely ill-conditioned and computationally prohibitive. Some recent works propose 
to use a simplified version of Navier-Stokes equations to circumvent this issue (see 
e.g. [ini [3H])- On the other hand, instead of relying the prior on the Navier-Stokes 
equations, spatial regularizers of u have been proposed to serve as a prior. A first 
approach in this direction is to assume a low-dimensional parametric form for u, see 
e.g. [H [13]. A second strategy is to choose a prior that introduces some spatial 
smoothness for u. Typically, the regularisation penalises in this case the norm of 
the gradient (or higher derivatives) of u, see e.g. 1351 [S5] . Finally, a third 

approach consists in the introduction of a self-similar constraint on u. Self-similarity 
is a well-known feature of turbulence, theoretically and experimentally attested, see 
e.g. ^36j. An attempt in this direction has been conducted in j221 [23]. Besides, a 
general family of self-similar regularizers has been introduced in [49]. In the same 
spirit, our choice for the prior of u is the divergence-free isotropic fractional Brownian 
motion (fBm), as we now justify. 

In addition to self-similarity, we assume u to be Gaussian. Non-Gaussian turbu- 
lent fields is an interesting alternative, which could potentially describe more accu- 
rately the structure of turbulence, but they will not been considered in this paper. 
Note that the definition of such models is still an active domain of investigation, see 
O [571 US 113 • Assuming further the stationarity of the increments of u leads nec- 
essarily to u being a vector fBm (see [48]). In order to satisfy the incompressibility 
constraint, u is finally demanded to be divergent-free. Although divergent-free fBm's 
can be viewed as a limiting case of the vector fBm's considered in [?5], we provide 
a proper definition in section ^ In particular, a spectral integral representation is 
deduced. 

This choice of prior involves in practice fractional Laplacian operators, which 
are numerically delicate to implement. Indeed, there is a lack of effective algorithms 
in the literature able to deal in practice with those particular priors. In |49| . the 
authors circumvent this issue in their numerical applications by limiting themselves 
to non-fractional settings. To tackle this problem, we propose to decompose fBms 
on well-chosen wavelet bases. This strategy allows us to expand fractional Laplacian 
operators on the wavelet components and makes their computation feasible. We focus 
on two particular wavelet bases. 

- As a first alternative, we propose to design wavelets as whitening filters for 
divergence-free isotropic fBms. For scalar fBms indexed by time, fractional 
wavelet bases represent ideal whitening filters [TT] [3S] . In [iSl [JT] [15] , the au- 
thors extend these bases to the case of isotropic vector fields, namely they use 
the so-called fractional Laplacian polyharmonic spline wavelets. In our case 
of divergence- free isotropic vector fBms, we show that for any mother wavelet, 
fractional Laplacian wavelet series composed with the Leray projector is an 
appropriate whitening filter. 

- The second alternative is to use divergence-free wavelet bases, which are well- 
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suited to our case. These wavelets simplify the decomposition, since they do 
not involve fractional operators and make superfluous the use of Leray pro- 
jector [131 [HI IMl US] ■ However, the wavelets coefficients are then correlated. 



As seen in section 5.1 implementation for the first bases can be accurately performed 
in the Fourier domain. For the second bases, the computation of the associated pos- 
terior density relies on the covariances of the wavelet coefficients. We provide a closed 
form expression for these covariances, which allows us to propose an approximation 
of the posterior using wavelet connection coefficients, see sections |5.2.1| and [5.2.2[ 

Moreover, an additional pleasant feature of wavelet representations is that it is 
well adapted to the non-convex optimization problem of MAP estimation, since it 
naturally provides a multi-resolution approach. As a matter of fact, this approach 
has proved to be experimentally efficient for motion estimation problems (subjected to 
different priors) 13, 25 , 50 . The optimization algorithms rely on Fast Fourier Trans- 
forms (FFT) or on Fast Wavelet Transforms (FWT). Finally, for the divergence- free 
wavelet decomposition, we introduce an approximated MAP optimization procedure, 
which turns out to be by far the fastest algorithm while being as accurate as the other 
approaches. 

A numerical comparison with state-of-the-art procedures is performed in section [6] 
on a general benchmark of divergence-free fBm. Note that, while these fields match 
perfectly our priors, they are likely to be in agreement with real fluid flows according 
to turbulence phenomenology, as mentioned before. As a result, our procedure seems 
to recover more accurately the hidden motion field, in particular when the Hurst 
parameter H corresponds to 2D and 3D turbulence models, respectively H — 1 and 
= I according to physics. 

The paper is organized as follows. In section [2] we give a spectral representation 
of divergence-free isotropic fBms, which will serve as the basic definition all along the 
paper. In section [3| we then provide the two wavelet representations described above. 
Section [4] displays the Bayesian modeling of turbulent flows in terms of wavelet coef- 
ficients and considers the associated MAP estimators. In section [5| gradient descent 
optimization algorithms achieving MAP estimation are presented for both wavelet 
bases. The numerical evaluation of the proposed regularizers is presented in sec- 
tion [6] Finally, the appendix gathers the technical proofs and some details about 
the algorithms and the computation of fractional Laplacian wavelet connection coef- 
ficients. 

2. Spectral definition of divergence-free isotropic fBm. The unique self- 
similar zero-mean Gaussian process with stationary increments is the fractional Brow- 
nian motion (fBm) , introduced in J34^ . The definition of scalar fBm indexed by time 
has been extended: to the case of multi-dimensional state spaces, i.e. vector processes 
indexed by time [5] ; to the field case, namely scalar isotropic fields [41' ; and to both 
cases, i.e. isotropic vector fields (551 148| . 

As far as we are interested in the construction of a prior for turbulent vector fields, 
we are particularly concerned with divergent- free isotropic vector fBm in the following 
sense. For sake of clarity, this work is restricted to the bi-variate case although it is 
possible to extend our results to higher dimensions. 

Definition 2.1. A bi-variate field u(x) = (ui(x), U2(x))^, x G M^, is a 
divergent-free isotropic bi-variate fBm with parameter < H < 1 if 

• u is Gaussian; 

• u is self-similar, i.e. for any A > 0, u(Ax) = A^u(x),- 
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• u has stationary increments, i.e. Vh G M , u(x) — u(x — h) is stationary; 

• u is isotropic, i.e. for any rotation matrix M , u(x) = u(Mx) ; 

• u is divergent-free, i.e. div u = almost surely. 

Since a fBm is almost surely not difFerentiable, the divergence operator in the last 
property is to be understood in a weak sense: for any test function ip e C^(IR^) with 
a fast decay at infinity, 



div u = (div u, ip) 



(u/W), 



where (., .) (resp. (./.)) denotes the inner product of two scalar (resp. two bi-variate) 
functions, and V denotes the gradient operator. 

Isotropic fractional Brownian vector fields have been introduced by Tafti and 
Unser |3S] as the weak solution of a fractional Poisson equation, or whitening equation: 



H + 1 



(2.1) 



where (— A)^ ^ is a fractional Laplacian operator depending on the Hurst parameter 
H > and some bi-variate complex constant ^ = (^1,^2), cr is a positive constant and 
is a vector of two independent Gaussian white noise. This gives rise to the solution 



a(-A) "^V, 



(2.2) 



where (— A)_^ ^ is the inverse of the operator involved in (2.1). The equation (2.11 



and the solution (2.2) are to be interpreted in the sense of generalized random fields. 



see [TB] and [H] , which explains the term "weak solution" . We refer to [17] and [15] 
for more details, including the definition of the above-mentioned operators. 

Isotropic divergence-free fractional Brownian vector fields, as we require, can be 
constructed as a limiting case of the solution (2.2 1, namely when ^1 = cx), see [48] 
section 4.7. It is worth mentioning that this case indeed corresponds to a limiting 
case of equation (|2.1|). In fact this equation has no sense for ^1 = 00, but the solution 



in (2.2) is well-defined when S^i — 00. 



In the following proposition, we provide a spectral representation of the isotropic 



divergence- free fractional Brownian vector fields when < < 1, as defined by (2.2 1 
when ^1 = 00 and ^2 = (the last condition is just a normalization so as to identify 
the constant a). This proposition gives an explicit representation, therefore it can 
also been considered as an alternative definition. As a matter of fact, we will only 
refer in the sequel to this representation when we consider isotropic divergence-free 
fractional Brownian vector fields. 



Proposition 2.2. The isotropic divergence- free fractional Brownian vector field 
u, as defined by ( [2.2[ ) when = 00 and ^2 = admits the following representation, 
for anyHe (0,1), 



u(x) 



a 
2tt 



(e 



1) 



IkF+i 



kk^ 



W{dk), X e 



(2.3) 



where I denotes the identity matrix and W = (W"i,W^2)"^ denotes a bi-variate stan- 
dard Gaussian spectral measure, i.e. Wi and W2 are independent and for i — 1,2, for 
any Borel sets A, B in Wi{A) is a centered complex Gaussian random variable, 
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W^{A) = Wi(-A) andE{W^{A)W^{B)) ^\AnB\. 

As a by-product, we deduce the following structure matrix function that charac- 
terizes the law of the Gaussian vector field u: for any i,j — 1,2, for any xi, X2, X3, X4 

in M2 



Sij(xi,X2,X3,X4) = E[(Ui(x2) - Ui(xi))(Uj(x4) - Uj(x3))] 

= a^H (i^f (X2 - X3) - (x2 - X4) - (xi - X3) + J^f (xi - X4)) , (2.4) 
where ch = r(l - H)/{tt2^"+^T{H + 1)H{2H + 2)) and 

^(2H 



F«(x) 



1)I-2H-. 



In particular, taking xi = X2 — h and X3 = X4 — h, for some h G M^, shows the power- 
law structure for the second order moment of the increments. We can also deduce 
the generalized power spectrum Ej, j = 1,2, of each component Uj of u, defined 
implicitly by (see [HI HI]): 



Sji(xi,X2,X3,X4) 

We deduce from the proof in Appendix |A.1[ 



ik-(x2-X3) _ gik-(xi-X4) _|_ gik-(xi-X3)\ 



E,(k)=a^ 1 



-2H-2 



(2.5) 



While all properties required in Definition 2.1 can be found in [58] going back 



to the definition (2.2) with ^1 =00, they are straightforward consequences of the 



spectral representation (2.3): Gaussianity and i7-self-similarity of u are easily seen 



from (2.3); Stationarity of the increments and isotropy follow from (2.4); Finally the 



divergence- free property of u is a consequence of the presence of the Leray projection 
operator I — -n^rj in (2.3), and will appear clearly in the wavelet decomposition of 



u considered in the next section. 



Remark 1. The definition ofu in (2.2 1 is not only valid for H € (0, 1) but also for 



H > 1, see 148^ - In- the same way, the spectral representation ( |2.3[ ) can he extended to 
H > \ by application of successive integrations as in 141^ , leading to the representation 



u(x) = 



2n 



CT f (e^''---ES(*k.x)Vj- 



Ik||^+i 



kk 



T 1 



IkiP 



VK(dk). 



(2.6) 



In this case u has stationary increments of order N = [77J +1, in the sense that 
for any (h^v, • . • , hi) e (M^)"'^, the successive symmetric differences D^^^ . . . _DhiU(x) 



form a stationary process, where : f(-) 
write specifically in this case: 



f(- + ^) ~ f(' ~ §)• These increments 



-Dh„ • . .£'hiu(x) 



a 
2^ 



N 



„ik-: 



I sm 



kh. 



kk^ 



W{dk). 
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3. Wavelet representations of divergence-free isotropic fBm. This sec- 
tion presents wavelet expansions of the divergence-free isotropic fBm representation of 
proposition [2?2] As pointed out in the introduction, wavelet bases are chosen in order 
to make fractional calculus involved by fBms feasible and effective. In practice we 
observe u on a bounded domain, say [0, 1]^. For this reason, we focus on the wavelet 
decomposition of u in (L^([0, 1]^))^. In the first part, the expansion of u relies on a 
fractional Laplacian wavelet basis and finally involves uncorrelated bi-variate random 
coefficients. In the second part, the expansion relies on a divergence-free wavelets 
basis, which is well-adapted to our setting, and involves mono-variate but correlated 
random coefficients. These two wavelet expansions will then allow us to derive ef- 
ficient optimization algorithms (see section [s]) solving the MAP estimation problem 
presented in the introduction. 

The construction of these two decompositions relies on an orthonormal wavelet 
basis of i^([0, 1]). There are several methods to build orthonormal wavelet bases of 
L^([0, 1]) (see chapter 7). For ease of presentation, we consider the simplest one, 
which consists in periodizing scalar wavelets of LP' (M) . Let be a mother wavelet with 
compact support and its associated wavelets dilated at scale 2^*+^ and translated by 

^,,,(a;)^2(^-i)/V(2^-^x-^). (3.1) 

The wavelet set {'4>(,^s{x)]x S M,^, s G Z} form an orthonormal basis of L^(M). The 
periodized wavelets V'^'T' ^, s G ^ is then defined by 

k —+oo 

C(^)= E ^^A^ + k), xe[0,l]. (3.2) 

k— — oo 

The set {ipYs ^ ^ [0, 1], s > 0, < £ < 2''^^} with the indicator function over 
[0, 1] form an orthonormal basis of i^([0, 1]). We will assume in the sequel that ip has 
M > H vanishing moments and is H + 2 times differentiable. These conditions will 
allow us to consider fractional derivatives and integrations up to order H + 2. 

3.1. Fractional Laplacian wavelet decomposition. An isotropic divergence- 
free fBm can in principle be decomposed in any wavelet basis of (L^(IR^))^. However, 
this kind of decomposition typically involves a sum depending on arbitrarily large 
scales, which is not suitable for application (see [3S] in the standard fBm case). In 
contrast, we show in Proposition |3 . 1 1 below that in our case where u is considered on 
the compact domain [0, 1]^, and under a mild condition (namely Jj^ u(x)(ix = 0), 
the divergence-free fBm enjoys a simple tractable decomposition in (i^([0, 1]^))^ with 
respect to fractional Laplacian wavelets. 

The bi-dimensional wavelet basis of (i^([0, 1]^))^ is constructed from i/;^'^^ as 
follows. Denoting the indicator function over the set A, we form the three following 
wavelet sets: 

{<i>£i,.iA0 = C'i,(a:i)I[04](^2); 0<4 <2^^-\ Si > 0}, 
{'i>ofi,i,,s,=l[o,i]{xiWCJx2); 0<£2<2''-\s2 > 0}, 

= C'.i(^i)Cr.2(^2); o<h<r'-',o<i2<r^-\si,s2>o}. 

Let us denote by il the set of indices {£,s) — (^i, si,^2, S2) involved in these three 
sets and {^e.s', (-^i s) G il} the union of them. An orthonormal basis of L'^([0, 1]^) is 
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finally the union of the latter set with the indicator function I[o^i]2(x) (see Theorem 
7.16 in [33] with J = 0). A basis of the product space (£^([0, i]^))^ is deduced by 
tensorial product. 

Let us now introduce the fractional Laplacian wavelets. They correspond to an 
integration of order iJ + 1 of ^g,s- 



(x)A 



1 



i2ny 



ik-x I 



-H-1 



B2 



(3.3) 



where $ denotes the Fourier transform of $. This fractional integration operator can 
be denoted in the spatial domain by (—A) s ^ following leading to the relation 



''(x) = Note that / {i, s) e n} constitutes a new 

family of wavelets, which is not orthogonal, unlike {^e.s; {i, s) £ ft}, but biorthogonal, 

where ^^'^'''^^ is the dual wavelet of i.e. {^["g^^^ -.^[vf,^^^) = S^^g'Ss.s' for 

all €,/,s',s' (see [I1IH1I35]). 

Let us finally recall the definition of the Leray projector, denoted by V, that maps 
square-integrable bi-variate functions v onto the space of divergence-free functions: 



i-H-l), 



1 



(2^)2 





■ kk^ " 


Is? 





v(k)dk. 



(3.4) 



This projection operator can be represented in the spatial domain by 7'(v) = [v — 
Vi(V.v)]. 

We are now in position to present the wavelet decomposition of u in (L^([0, 1]^))^- 
We assume for simplicity that /jp u(x)(ix = 0. Note that in practice, we observe u 
on a lattice with n x n sites, so the latter assumption roughly means that the mean 
value of u on this lattice is assumed to be negligible. 



Proposition 3.1. Let u be an isotropic divergence-free fBm with parameter 
H G (0, 1). Assuming J^^^ u(x)(ix — 0, we have for all x G [0, 1]^; 



u(x)= J2 ^ 
(e.s)en 



(-H-1) 



(x), 



(3.5) 



where coefficients e^^g o,re i.i.d. bi-variate zero-mean Gaussian random variables with 
variance (27r(T)^, and where ^gg^ ^\ (^i *) G ^^^6 defined in (3.3). 



Remark 2. As shown from the proof the simple form (3.5) holds because the 
orthonormal basis o/L^([0, 1]^) that we have considered (the periodic wavelet basis) 
involves a unique scaling function which is the indicator function I[o,i]2(x). It is there- 
fore clear that the representation (3.5) remains valid for any other wavelet basis of 
L^([0, 1]^) having the indicator function as its unique scaling function, 
particular the case for the folded wavelet basis 



This is in 



Remark 3. Following remark [7} it is easy to adapt the proof of Proposition \3.1\ 
to the case H > 1 and deduce that ( |3.5[ ) remains also valid for H > 1, since the 
assumption J^^ u(x)(ix ~ removes all constant terms in the wavelet expansion. 
Moreover, when H is an integer, there exists no clear representation of the divergence- 



free fBm. Since the representation ( 3.5 ) is well defined for any H , we choose to extend 
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by continuity (3.5) to integer values of H. 



3.2. Divergence-free wavelet decomposition. Since u is continuous and 

divergent-free, tlren u e L'^iy{[0, 1]^), where ^^^^([O, l]'^) denotes ttie space of divergence- 
free bi-variate fields in [0, 1]^, i.e. 

iL([0,l]')-{ve (L2([0,lp)f:divv = 0}. 

In ttiis section we propose a decomposition of u onto a biorthogonal wavelets basis 
of L'^iyilO, 1]^)- Such basis is constructed from an orthonormal basis of ^^([0, 1]^), as 
described below. 

Let us start from the periodized wavelets i,s G Z defined in (3.2). The 



primal divergence-free wavelets — i^l s^^e s)^ biorthogonal wavelets 



basis of L'^iyiiO, 1]^) are defined for {£, s) e by 



- for 0<£i<2^i-i,0 <£2<2''^-\si,S2 > : 

*t(x)=c:..(^i)c;(^2) and ^u^) ^ 

- for 0<4< 2^^1^1,51 >0, ^2=0,.S2 = 0: 

^IJ^) = and ^U^) = -Vn,(a:i)I[o,i](^2), 

- for 4 = 0,si = 0, 0<4<2^2-\s2 > : 

*L(x) -I[oa](2^i)Cr.(^2) and ^^ .(x) - 0, 

where Tp' denotes the derivative ofip. To complete the primal wavelets family, the func- 
tion 4'o(x) = (I[o,i]2(x),I[o_i]2(x))^ is superimposed, leading to the primal divergence- 
free wavelets family {"^e^s', {(■: •§) € riUO}. Note that we have = curl[$£ g] where 
curl = ( —-rp-Y and „ is the orthonormal basis constructed in section 

^ CtX2 ' OXi f 



3.1 



The dual wavelets ^ — ^\ s-'^\ s)^ ^'^^ constructed in order to be 

biorthogonal to \I'£,s, i-e- {'^ i.s / e .s') = S^^^'Ss^s' for all £,£',s,s'. They are given 
by *o = *o and for {£, s) e Q., 

- for 0<£i<2"i-i,0 <£2<2''^~\si,S2 > : 

^t(x)=-c:.,(xi) rr,::jx)dx and ^u^)=rCsM) r^zs^)dx, 

JQ JO 

- for 0<4<2"i-\si > 0, ^2 = 0,S2 = : 

§L(x)=0 and ^'2,(x)=I[o.i](x2) / 'v'rsi(^)^^' 

Jo 

- for £i = 0,si = 0, 0<^2<2'"^-\s2 > : 

*i,,(x) = -I[o.i](xi) / riZ2i^)dx and *^,,(x) = 0. 



The decomposition of u in the divergence-free biorthogonal wavelet basis de- 
scribed above writes, for all x e [0, 1]^, 

U(x) = ^ (i^,s*£,s(x), 

{e,s)enuo 
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where d^ s — {u/^e,s) are the random divergence-free wavelet coefficients. Unhke de- 
composition (3.5), these wavelets coefficients are in general correlated. The following 
proposition describes their covariance structure under the assumption Jjjj u(x)dx = 

0. Note that in this case, since (u/^o) = 0, the decomposition of u becomes 

u(x)= J2 dt.s'^iA^)- (3.6) 

Proposition 3.2. Let u be an isotropic divergence- free fBm with parameter 
H G (0,1). Assume that Jjp u(x)(ix — 0. Then, the coefficients di^s, (-^jS) e fl 
are zero-mean correlated Gaussian random variables, characterized by the following 
covariance function: for any {£,s), {£',s') E ft, 

S(€,s,€',s') ^ EK,,d,,,,,] = (27ra)2($(/-2),ci>(7^«-2)), (3.7) 

where "^^ is defined in (3.3). 

In practice, the fBm u is represented by a finite number of wavelets coefficients 
dg s, so that the covariance function s, s') becomes a matrix, say S„ (see 
section 4.2). The following proposition will allow us to construct an approximated 



inverse matrix of S„ for later use. 

Proposition 3.3. The operator's : £'^{fl) £^{n) defined for any a £ £'^{n) by: 

[Sak. - ae',s''^{£,s,£'.,s'), V(£, s) e f!, (3.8) 

{t,s')en 

admits an inverse operator given for any a G £^[n) by: 

[S''ak«= (^t,s'^'^{£.s.,£',s'), yi£,s)en, (3.9) 

{e',s')en 

where for any {£, s), {£' , s') G f2, 



In the paper remainder, we will assume Jj^ u(x)dx = to fit the assumptions 
of Propositions |3.1| and |3.2[ 

4. MAP estimation. At this point, we have gathered all ingredients necessary 
to formalize properly the problem of estimating the incompressible velocity field u of 
the introduction by a MAP approach. In this section, we ffi-st justify the choice for 
the likelihood model psy\u from physically sound consideration. Then, we express the 
fBm prior for u in terms of the two wavelet expansions presented before. Finally, we 
deduce the optimization problems to solve in practice in order to obtain the MAP 
estimates. 
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4.1. Likelihood model. Solving a turbulent motion inverse problem consists 
in recovering a deformation field, denoted later on by u, from the observation of an 
image pair. To this aim, the estimation classically relies on a likelihood model linking 
the unknown deformation to the observed image pair {yoij/i}- We first assume that 
this image pair is the solution y(x, t) e C^([0, 1]^ x M) taken at t = and t = 1 of the 
following transport equation, so-called in the image processing literature "optic-flow 
equation" : 

^(x,t)-fv(x,i)-Vx2/(x,i)=0 ^^^^ 
y(x,0) = yo(x) 

where v(x, i) is the transportation field that verifies, at any time f G M, v(.,i) g 
L%y{[0, 1]^) due to the incompressibility constraint. It is well known [JD] that we can 
write ?/o(x) = ?/i(Xj(x)) where the function t Xg(x), known as the characteristic 



curves of the partial differential equation (4.1), is the solution of the system: 

|x*(x)=v(X^(x),i) 
XS](x)=x 



Let us remark that system (4.1) constitutes a relevant model frequently encountered 
in physics, in particular in fluid mechanics: it describes the non-diffusive advection of 
a passive scalar by the flow v [5D]. Assuming that v(Xq(x), s)ds = v(x,0), then 
denoting u(x) = v(x, 0) we obtain by time integration the so-called Displaced Frame 
Difference (DFD) constraint: 

yo(x) = yi(x-|-u(x)). 

However, in the context of a physical experiment, the observed image pair is likely 
to be corrupted by noise measurement errors. In this context, a standard approach 
to estimate u is to adopt a probabilistic framework. We make the assumption that 
the quantity j/i(x 4^ u(x)) — ya(x) is corrupted by a centered Gaussian noise. More 
precisely, we assume that the so-called data-term DFD functional 

5y(u)^^|lyi(x + u(x))-yo(x)|P; u€LL([0,in, yo, yi e ^'([0, 1]^), (4.2) 
follows an exponential law, so that the likelihood model writes 

P5,|u = /?exp-^*^W, (4.3) 

where /3 is a positive constant. 

A straightforward criterion to estimate deformation field u from the observation 
of {yo,yi} is to maximize the likelihood: argminQg^2 Sy{u). Although searching 
the deformation field u in L^^^ instead of (L^)^ reduces by a factor two the number of 
degrees of freedom, this problem remains ill-conditioned. This is due to the so-called 
aperture problem |25l I52j . Therefore, a Bayesian framework introducing prior infor- 
mation for u is needed for regularization of the solution. 
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4.2. Prior models. As explained in introduction, we choose as a prior for u the 
divergence-free fBm described in section [2j From the wavelets decomposition (3.51 or 
(3.6), an isotropic divergence-free fBm prior can 

(i) either be represented by a Leray projection of a fractional Laplacian wavelet 
series, whose coefficients are distributed according to independent standard 
normal distributions. 

(ii) or by divergence-free wavelet series, whose coefficients are correlated accord- 
ing to (3.7). 

In practice, the images j/o and yi are of size nxn pixels. Accordingly we will truncate 
the series (3.5) and (3.6) with respect to the index set 

rin ^ {i,s e il; si, S2 < Sn = log2in)}. 

In other words, wavelet coefficients revealing scales smaller than the pixel size will 
be neglected. We thus consider the vector of fractional Laplacian wavelet coefficients 
e„ = (e^, e^)^ with e'^ = {e^ (£, s) e fi„} and the vector of divergence-free wavelet 
coefficients d„ = {d^,^; {f;^) & ^n}- 

In order to express the prior for u in this finite dimensional case, we need to 
introduce some notation. Let us first denote by S„ the matrix whose element at row 
{£,s) and column {£\s') is E[(i^_s'^f',s']j that is the covariance matrix of c?„. From 
proposition 3.2 we have for all {£,s), {£' , s') S ^In- 



i-H-2} ^{-H-2) 



Its inverse may be approached by the matrix denoted by S„ ^ composed of elements 



(27rCT)2 



1 



(4.4) 



From proposition 3.3 this approximation becomes accurate for n sufficiently large. 
Second, for any H > Owe introduce operator ^i^^^^^ from (^^(ri„))^ into (i^([0, 1]^))^ 
defined for any e„ e (£^(il„))^ by 



(^,s)en„ 



(4.5) 



We introduce analogously operator ^„ : ^^(ri„) — >■ ^^^^([0, 1]^) defined for any d G 
eH^n) by 



(f,s)en„ 

We are now in position to state the two representations of the prior for u: 
(i) based on fractional Laplacian wavelet decomposition (3.5): 



(4.6) 



u = vm 



i-H-l] 



e„), with pe„ = 



1 



(27r)"/2(27 



_g 2(2ircr)'i 



(4.7) 



(ii) based on divergence-free decomposition (3.6): 

1 



with — 



(27r)t det5(S„ 



(4.8) 
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4.3. Maximum a posteriori estimation. From (4.7 1 or (4.8), u reduces to 



the knowledge of wavelet coefficients e„ or d„. Let us rewrite the likelihood in terms 
of these coefficients: 



P5y\e„ = PSy\d^ = P exp 



-0Sy 



with the DFD data term (4.2 1 rewritten as 



5y{-) = ^l|yi(x, •) -2/o(x)|P, 



(4.9) 



where 



yi(x,e„)^yi (x + 7'(*i-^-i)e„)(x)) or yi(x, d„) ^ (x + *„d„(x)) . (4.10) 



The MAP estimates are defined by: 



e„* =argmaxp5j^|g^Pe„ and d*^ ^ &Ygu\&yipsy\d„Pd^- 



(4.11) 



Solving the MAP problems (4.111 is equivalent to minimize minus the logarithm of 
the posterior distributions: 



i) with respect to fractional Laplacian wavelet coefficients 

e„* = argmin{(52/(e„) +7^(e„)} 
1 



7e(e„) ^ 



T 



2/3(27r(T)2 " 

ii) with respect to divergence-free wavelet coefficients 



d* = aigim\i{5y{dn) + 7^(d„)} 

1 ' 



Ti-idn) = -^dn S„ dn, 



(4.12) 
(4.13) 

(4.14) 
(4.15) 



where TVs are so-called regularizers. 

In the following we will assume that the Hurst exponent H is known. Posterior 
models described above are thus characterized by two free-parameters, namely j3 
and a. However, MAP estimates (4.12) and (4.14) only depend on the product of 

This 



these two parameters, forming the so-called regularization parameter 



regularization parameter explicitly appears in (4.13) while it is partially hidden in 



(27I-0-)2 



the covariance matrix inverse in (4.15). Let us mention that an empirical study 
shows a low sensitivity to the choice of the regularization parameter (see section [6]). 
Nevertheless, estimation techniques exist when no prior knowledge is available for the 
adjustment of Hurst exponent or regularization parameter [2 H l32 j |44]. 



5. Optimization. In this section, we introduce algorithms to solve (4.12) and 



(4.14). To deal with these non-convex optimization problems, we use a Limited- 
memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) procedure, i.e. a quasi-Newton 
method with a line search strategy, subject to the strong Wolf conditions [37]. More- 
over, as suggested in |13j we propose to enhance optimization by solving a sequence of 
nested problems, in which the MAP solution is sequentially sought within higher reso- 
lution spaces. More precisely, wavelet coefficients are estimated sequentially from the 
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coarsest scale 2" to the finest one 2 At each scale, problems (4.12) and (4.14) are 



solved by the LBFGS method with respect to a growing subset of wavelet coefhcients. 
This subset includes all coefficients from the coarsest scale to the current one 2"**, 
coefficients estimated at previous coarser scales being used as the initialization point 
of the gradient descent. This strategy enables to update those coarser coefficients 
while estimating new details at current scale 2~*. 

To implement the above procedure, the functional gradient and the functional 
itself in (4.12) and (4.14) need to be evaluated at any point e„ and d„. Note that 



once the functional gradient is determined, the functional value is simple to deduce: 



the evaluation of (4.10) needed by (4.9) will be a precondition to the computation of 



the DFD functional gradient while (4.13 ) or (4.15 ) may be derived from their gradients 
by simple scalar product with the vector of wavelet coefficients e„ or (i„. 

The next section provides algorithms to compute exactly the functional gradient 
for the two different wavelet representations. Besides, an approached gradient com- 
putation is proposed to accelerate the algorithm in the case of the divergence-free 
wavelet representation. 

5.1. Projected fractional wavelet series. 

Proposition 5.1. Let Vi'v denote the i-th component of the projection of a 



bi-variate vector -v onto L~^^^{[0^1\ ). The gradient of functional minimized in (4.12) 
with respect to e\ ^ is d^i^ Sy{en) + d^i^ Ti-i^n) where for i = 1, 2 and for all (£, s) e ri„; 



d,.Jy{e,,) = ((-A)"^7'.[(yi(-,en)-yo(-))Vyi(-,e„)](x),<i>,,,(x)) (5.1) 
1 



(5.2) 



Based on these formula, we derive a spectral method for the computation of 
the gradient of functional minimized in (4.12). It is based on FFT and FWT with 
recursive filter banks: 



Algorithm 1. (functional gradient w.r.t e„) 

i) compute the FFT of the components of (j/i — j/o)Vyi, 

ii) compute fractional differentiation and Leray projection in Fourier domain 
Hi) compute the inverse FFT to get {~A)^^Vi[{yi{-,en) — yo{-))'Vyi{-,£n)] 
iv) decompose each component by FWT using the orthogonal wavelets s to 
obtain the data-term gradient (5.1), 



v) Derive functional gradient by adding vector (5.2) to the data-term gradient 



In order to evaluate t/i or its gradient in the above algorithm, one needs to reconstruct 



the unknown u (4.7) appearing in (4.10) from the fractional Laplacian wavelet coef- 



ficients e„. This can be done by the following spectral computatior[^ of the inverse 
fractional wavelet transform. Indeed, by commuting Leray projector with fractional 
integration, (4.7) can be rewritten as: 



u = 



1 



dke' 



— ik-: 



(0). 



(5.3) 



^Let us remark that this reeonstruetion essentially differs from a direct spectral fBm generation 
method which is known to bring aliasing side-effects, see [3]. 
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where we recall that ^^"^ is defined in (4.5) 



This yields the following reconstruction algorithm: 

Algorithm 2. (reconstruction of fBm from e„) 

i) reconstruct ^n'^en from e„ by inverse FWT of each component using orthog- 
onal wavelets {£, s) e iln}, 
ii) compute FFT of the two components of the latter function, 
Hi) compute Leray projection and fractional differentiation in Fourier domain, 
iv) compute inverse FFT of the two components to obtain u. 

Algorithms [T] and [2] yield the ingredients necessary to approach the MAP estimate 
e* with a gradient descent method of theoretical complexity bounded by the FFT 
algorithm in 0{n log n). Nevertheless, the computation bottleneck comes mainly from 
the number of transforms that are required at each gradient decent step: 4 FFT, 4 
inverse FFT, 2 FWT and 2 inverse FWT. 

5.2. Divergence-free wavelet series. 

5.2.1. Exact method. For notational convenience we define operators ^J^^^'^^^-' 
€^(n„) — >• i^([0, 1]^) for i = 1,2 as the two components of the operator defined in 



(4.5): 



^^l-"-'^e,,. (5.4) 



1 

—n n 
$2,(-H-l) 2 
—n n 

Note that we have = 

Proposition 5.2. The gradient of functional minimized in (4.141 with respect 
to dg^s is dd^ Jy{dn) + dd^^Tl{dn) where for all {£, s) e ^n-' 

dd,Jy{dn) = {{yi ~ z/o)Vyi/*^,,) (5.5) 

and 

dd.,.ndn) = Jli^(^n^"^'^d^^ $2+'^ (5.6) 
Based on the previous result, we derive a spectral method for the computation of 



the gradient of functional minimized in (4. 14 1. It is based again on FFT and FWT 
with recursive filter banks: 

Algorithm 3. (functional gradient w.r.t d„) 

i) decompose (yi— 2/o)Vyi by FWT using dual divergence-free wavelets {£, s) € 



fin} to obtain the data-term gradient (5.51, 
ii) compute inverse FWT of dn using orthogonal wavelets {$£.s; (-^j € ^n}- 
Hi) compute FFT and apply operator A^^^ in Fourier domaii^ 
iv) compute inverse FFT to get A^^'^ {^l^^^^ dn) ■ 



v) compute FWT using orthogonal wavelets {^£,s'i (-^j *) € ^n} get (5.6) 



vi) derive functional gradient by adding (5.6) to the data-term gradient. 



^ This supposes 2H+A differentiability of ip. If we only assume H+2 differentiability, algoritlim|3] 
can be modified as described in appendix IbI 
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Reconstruction of u from coefficients dn is needed to compute yi . It is easily done 
by an inverse FWT. 

Algorithm 4. (reconstruction of fBm from d„) 

i) Reconstruct u = ^^dn from dn by inverse FWT using divergence-free wavelets 

Algorithms [3] and |4] yield the ingredients necessary to approach the MAP estimate 
d* with a gradient descent method of theoretical complexity bounded again by the 
FFT algorithm in 0{n\ogn). However, it is less time-consuming than the approach 
based on fractional Laplacian wavelet series since it requires for each gradient decent 
step: 1 FFT, 1 inverse FFT, 2 FWT and 2 inverse FWT. Moreover let us remark that 
a simple FWT with recursive filter banks corresponding to the dual divergence-free 
wavelet basis |29j is required to compute the data-term gradient. 

5.2.2. Approach method. In this section, we approach the divergence-free 
regularizer gradient (5.6) in terms of matrix products, which will make computation 
of the gradient more efficient since we avoid intensive use of FFT and FWT as in 
algorithm [sj From proposition 5.2 and the Parseval formula, we have for all (€, s) € 



^ (k),iikir+2$,,,(k)) (5.7) 



where ^ denotes the two-dimensional Fourier transform of ^g^s- We hereafter derive 
a separable approximation of (|5.7|) in terms of one-dimensional connection coefficients 



fe^As defined for all {£', s'), {£, s) e as 

V^f/) for < s,s' < s„, < £<2'*-\0 < f <2'*'-^ 

y. 

for s = 0,0 < s' < s„, l = Q,0< /<2"'-i 



Jl'.s'. 



(I[0,1], 

(I[o.i], 



92 



-92 



V^rj) for < s < s„, s' = 0, < ^ < 2 



'-^ r = 



^[0,1] 



for s = 0, s' = 0, Z = 0, = 



(5.8) 

Note that for any fixed a < (resp. a > 0), fi?\i^g exists whenever -0 possesses 
sufficient vanishing moments (resp. is sufficiently differentiable) . 

The two-dimensional scalar products in (5.7) can be written ^'(k), ||k||^''^^^''$^.s(k)). 
If G N, Newton's binomial theorem applies: 



\2{H+2) 



where 



H + 2 

i 



= {H + 2){H + 2- 1)...{H + 2-1 + 
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So if iJ e N, plugging (5.9 1 into (5.7) shows that ddi^'R{dn) can be expressed in a 
separable form: 



1 



H+2 



H + 2 



EAH+2-i) , f(i) 



12 ■S2 



(5.10) 



The latter formula can be expressed in terms of matrix products. Let F*^"^ be 
the matrix of size n x n composed at row index (£' , s') and column index (£, s) by the 
element /^^"]/ £ g. Let the n x n matrix whose element at row index {£i, si) and 
column index (^2iS2) is d^^s- Equation (5.10) can then be written 



1 



'H+2 



(5.11) 



2/3(27rcr)2 

where [MJ^.g denotes the (£, s)-th element of matrix M. 

In the general case where H ^ N, (5.11 ) does not hold anymore, but in the same 
spirit we consider the following natural approximatiorj^ 



ddi„'R-id„) 



2/3(27rcr) 



\H+2i 

E 

i=0 



H + 2 
i 



f(«+2-^)"wjfW + fW' 



(5.12 



where [H + 2J denotes the integer part of i7 + 2. 

We approximate matrices F*-"' by an off-line FWT of scaling function connection 
coefficients, as explained in appendix [C] and summarized in the following algorithm. 
These scaling function connection coefficients are in turn easily computable as a so- 
lution of a linear system, by adapting Beylkin's fast methods [71[S1|5], see appendix 



Algorithm 5. (ofT-line computation of matrices F^")) 
i) Compute scaling function connection coefficients 



i2 \ " 



^{2'"x-£')}t 



for any integers £,£' by inversion of linear system ( C.13[ ) (see the explicit 
solution (C.15)j, where ip denotes the scaling function associated to wavelet 
^. 

ii) Construct the discrete (2'^")-periodic function defined for any integers £,£' 
with 0<£,£' < 2"" by: 

for < |f -^1 < 2^"-i, 
[ei"^(£,f -2^") for 2""-'^ < \£' - £\ < 2'" . 



^This approximation may be justified by considering a truncated version of the Newton's gener- 
alized formula. However, a control of the approximation error of (??) is a technical issue, which is 
out of the scope of this paper. 
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Hi) ApproachF^"'^ by performing the (discrete) FWT of the latter two-dimensional 
function multiplied by factor 2*^" . 

Besides, note that for matrices F^") with a G N, connection coefficients computed 
in step i) of algorithm [s] are classicaUy obtained by solving an eigenvalue problem, see 



(4.14) 



p.2]. Using (5.12), we hereafter propose a fast algorithm for the minimization problem 



Algorithm 6. (approach functional gradient w.r.t c?„) 
i) Decompose (yi — yo)yyi by FWT using the dual divergence-free wavelet basis 



{^£,s;(-^iS) G ^n} to obtain the data-term gradient (5.5), 
ii) Derive functional gradient by adding the first terms of (5.; 
to the truncation chosen above) to the data-term gradient. 

The reconstruction algorithm is identical to algorithm [4j 

Algorithm [6] can be substituted to algorithm |3] to approach the MAP estimate d*. 
In theory, it requires a theoretical complexity bounded by matrix multiplication, i.e. 
0{n^). In order to reduce this theoretical complexity, one may take advantage of 
the very sparse structure of matrices F'"-*, which could be investigated as in |39] . 
However, in practice FFT and FWT are the bottleneck of the computational cost. 
Algorithm |6] requires only one FWT for each gradient step, in contrast to the numer- 
ous FFT and FWT involved in algorithm |3] All in all algorithm |6] turns out to be 
much faster than algorithm [3j 



6. Numerical evaluation. In this section, we assess the performance of the 
proposed estimation algorithms and compare them to recent state-of-the art methods. 
As explained below, the numerical evaluation relies on a synthetic benchmark of noisy 
image couples revealing the turbulence transport of a passive scalar. Turbulence is 
mimicked by synthesizing fBms for various Hurst exponents, which include H = ^ 
and H = 1 modeling respectively 3D and 2D turbulence. 

6.1. Divergence-free isotropic fBm fields generation. Divergence-free isotropic 
fBms were generated by a wavelet-based method relying on reconstruction formula 



(5.3). More precisely, in agreement with the fBm model ( |4.7[ ), the wavelets coeffi- 
cients {e^.s} were sampled according to standard Gaussian white noise. The fBms 
realizations were then synthesized by application of algorithm [2] 

In order to form an evaluation benchmark for the regularization model, a set of 
5 fBm realizations were synthesized (from an identical white noise realization) with 
a resolution 2~^ x 2~^ on the domain [0, 1]^. The wavelet generator was constructed 
from periodized Coiflets with 10 vanishing moments, so-called Coiflets-10 [33] • Note 
that these wavelets are 3 times differentiable so that our regularity assumptions are 
satisfied for any H < 1. The realizations are associated to Hurst exponents Hi — 0.01, 
H2 = H3 = ^, = I and = 1. Let us remark that the cases H = ^ and 
H = 1 are consistent with isotropic turbulence models introduced by Kolmogorov for 
3D flows |2Z| (resp. by Kraichnan for 2D flows I28j), while the case H = ^ constitutes 
a standard Brownian motion. Figure C.l displays vector fields u(x) — (ui(x), M2(x))* 



and vorticity maps 9j;U2(x) — dyUi{x.) associated to the 5 fBms. 
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Let the radial power spectrum of u be defined by: 

E{k) ^ I Si(x) + S2(x)dx, (6.1) 



where 5^ is the circle of radius n. It is easy to check that according to (2.5), this 
function is E{k) — ck~^^~^, where c > 0. Therefore, the spectra of our 5 different 
fBms decay exponentially with power respectively equal to —1.02, — |, —2, — | and 
-3. 

6.2. Synthetic image couple generation. To simulate the couple {yo,yi) in 



the data-term (4.2), we start by a fixed image yi and derive yo from the relation 



2/0 (x) = yi(x + u(x)) derived in section 4.1 Since x + u(x) does not necessarily 



lie on the pixel grid, we used cubic B-splines for interpolation of yi. To simulate 
realistic measurement conditions, the so-generated images were then corrupted by 
i.i.d. Gaussian noise yielding a peak signal to noise ratio (PSNR) on yo (resp. t/i) 



of 33.2 dB (resp. 33.5 dB). The resulting image pairs are displayed in figure C.2 for 



= I and iJ = 1. 

6.3. Optic-flow evaluation procedure. The divergence-free fBm fields were 



estimated according to a MAP criterion, solving minimization problems (4.12) or 



(4.14). The proposed approaches are compared to two other standard regularizers, 
which all require the choice of a basis to decompose u. In order to make relevant 
comparisons, we chose the divergence-free wavelet basis for all these alternatives. 
The wavelet generator was constructed from divergence-free biorthogonal Coiflets-10 
with periodic boundary conditions. Moreover, the optimization procedure used for all 
regularizers was the same and relied on an identical data-term. 

The five different estimation methods used for evaluation are listed and detailed 
hereafter. They are denoted A, B, C, D and E. Methods A and B are state-of-the-art 
algorithms while methods C, D and E implement the fBm prior. 

A - Penalization of norm of velocity components gradients. The most common 
approach in optic-fiow estimation, as first proposed in [24], is to penalize the 
norm of the velocity gradients. We used the wavelet-based implementation 
proposed in 

B - Penalization of norm of vorticity gradient. In fiuid motion estimation a 
popular approach is to penalize the norm of the vorticity gradient |lll 1451 
I52j . Here again, we used the wavelet-based implementation proposed in |25] . 
C - fBm regularization in a fractional Laplacian wavelet basis. This corresponds 
to solve (4.12) using algorithms [l] and [2] 



D - fBm regularization in a divergence-free wavelet basis. This corresponds to 
solve (4.14) computed without approximation using algorithms |3] and |4] 



E - Approached fBm regularization in a divergence-free wavelet basis. This cor- 



responds to solve (4.14) where the regularization term is approached using 
algorithms |4] and |6 

Each of the regularizer in methods A, B, C, D and E were optimally tuned, that is 
to say regularization coefficient were chosen (using a brute-force approach) in order 
to the RMSE detailed hereafter. Note that an implicit regularization by polynomial 
approximation has also been tested. It is a well-known approach in computer vision 
[11 [131 |3TJ [50j . The performances were clearly below the previous approaches, so we 
do not display the results in this paper. 

Let S denote the set of pixel sites. The two following criteria were used to evaluate 
the accuracy of estimated fields denoted by u*: the Root Mean Squared end-point 
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Error (RMSE) in pixel 



RMSE = Mr 7 |u* W - u(x; 



|2 



,2 ^ 

xe5 



and the Mean Barron Angular Error (MBAE) in degrees 



MBAE = ^ > arcos 



,2 ' 



^ V |u(x 



u*(x) • u(x) 

2 



where u represents the synthesized fBm. Moreover, we introduce a criterion to eval- 
uate the accuracy of reconstruction of the power-law decay of the radial power spec- 



trum. More precisely, in logarithmic coordinates the power-law (6.1) writes as an 
affine function of log(K) of the form — {2H + 1) log(K) + 7, depending on two parame- 
ters, namely the Hurst exponent H and the intercept 7, that can be explicitly related 



to a in (2.5 1. Performing a linear regression (by the ordinary least squares method) 
on the estimated spectrum in logarithmic coordinates, we obtain an estimation of 
the afhne function parameters denoted by H* and 7*. The quality criterion is then 
chosen to be the distance between the estimated and true afRne functions within 
the interval [log(Kmm), log(Kmaa;)] , called the Spectrum Absolute Error (SAE): 



(.log(K„aa;) 

SAE= / \2x{H - H*) +-f* -j\dx, 

Jlog(K„i„) 



In order to evaluate the power-law reconstruction at small scales, we chose K„im = 

10 and Krnax = n. 

Finally, we performed an additional visual comparison of the accuracy of resti- 
tuted vorticity maps. 



6.4. Results. In table of figure C.3 the performance of the proposed methods 



(C, D and E) can be compared in terms of RMSE, MBAE and SAE to state-of-the- 
art approaches (A and B). Let us comment these results. Methods C, D and E yield 
the best results with respect to each criterion. Method C, i.e. the method based 
on fractional Laplacian wavelets, provides the lowest RMSE and MBAE. An average 
RMSE gain of 19% with respect to the best state-of-the-art method is observed, with 
a peak at 26% for = i. However, considering the 3 criteria jointly, methods D 
and E, i.e. exact and approached method based on divergence-free wavelets, provide 
the best compromise. In particular, according to SAE it can be noticed that, unlike 
method C, methods D and E achieve to accurately reconstruct the power-law decay 



of the fBm spectrum. This is illustrated in figure C.4 Moreover, the approximation 
used to derive method E seems to be accurate since performance of E are very close to 
those of D. In figure [C?5| one can visualize estimated vorticity maps with the different 
methods for = | and H ^ 1, i.e. fBms modeling respectively 3D or 2D turbulence. 
This figure clearly shows the superiority of methods D and E in reconstructing the 
fractal structure of the vorticity fields. 



Plots of figure C.6 show the influence of the regularization parameter in terms of 
RMSE, MBAE and SAE for 7? = | and iJ = 1. Clear minima of RMSE and MBAE 
are visible for methods A, B, D and E. On the other hand, method C, i.e. fractional 
Laplacian wavelet basis, seems to be 'unstable' in the sense that it yields inhomoge- 
neous performances for small variations of regularization parameter. The saturation 
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of the RMSE and the MBAE for large values of the regularization parameter shows 
that, on the contrary to state-of-the-art methods, sensitivity of method D and E to 
the latter parameter is weak, in the sense that it yields reasonable estimation error 
even for regularization parameter far from an optimal value. This error saturation 



effect is illustrated in figure C.7 It displays vorticity maps produced by the different 



methods for a very large regularization coefficient. 

7. Conclusion. This work addresses the inverse-problem of estimating a hidden 
turbulent motion field from the observation of a pair of images. We adopt a Bayesian 
framework where we propose a family of divergence-free, isotropic, self-similar priors 
for this hidden field. Self-similarity and divergence-free are well known features of 
incompressible turbulence in statistical fluid mechanics. Our priors are bivariate frac- 
tional Brownian fields, resulting from the extra assumptions that the hidden field is 
Gaussian and has stationary increments. The main purpose of this article is the design 
of effective and efficient algorithms to achieve MAP estimation, by expanding these 
specific priors into well-chosen bases. From a spectral integral representation proved 
in proposition [2]2] we represent divergent-free fBms in two specific wavelet bases. The 



first option (proposition 3.1) is a fractional Laplacian wavelet basis which plays the 
role of a whitening filter in the sense that the wavelet coefficients are uncorrelated. 
The second alternative is to use a divergence-free wavelet basis, which is well-suited 
to our case. The latter wavelets simplify the decomposition, since they neither in- 
volve fractional operators nor Leray projector on the divergence-free functional space. 
However, the wavelets coefficients are then correlated. We provide a closed-form ex- 



pression for the induced correlation structure (proposition 3.2 1, which is necessary to 
implement this second approach in practice. For these two approaches, the algorithms 
to reach the MAP involve gradient based LBFGS optimization procedures and rely 
on fast transforms (FFT or/and FWT). Moreover we propose an approximation of 
the correlation structure of the coefficients in the divergence-free wavelets expansion. 
It is based on off-line computation of fractional Laplacian wavelet connection coeffi- 
cients. This approximation leads to the fastest algorithm without loss of accuracy. 
According to an intensive numerical evaluation carried out in section [6j all proposed 
algorithms clearly outperform the state-of-the-art methods. Finally, in the light of our 
experiments, the divergence- free wavelet expansion seems to be the most appropriated 
representation to solve our MAP inverse-problem. 
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Appendix A. Proofs. 



A.l. Proof of Proposition 2.2, Let us recall (see e.g. [5T]) that given a 
standard Gaussian spectral measure Wi, the integral / f{'k)Wi{dk) is well-defined 
whenever / G i^(M^), has zero expectation and for f,g in L^(IR^): 



E (^J f{k)W,{dk) J g{k)Wi{dk)^ ^ J f{k)g{k)dk. (A.I) 



In (Esl), the matrix P(k) = 



IlkiP 



corresponds to the Leray projection matrix 
in the Fourier domain. It is easily verified that all entries of P are in [0, 1]. For this 
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reason the integral (2.3 1 is well defined, since for all x e and all H E (0, 1), the 
function k i-^ (e*''''' - l)\\k.\\-"^'^ belongs to L'^iR"^). 



Let us show that the structure function of u is given by (2.4 1. For j — 1,2, 
denote the bivariate vector whose j-th component is equal to one while the other 



component is zero. For all i,j ~ 1,2, from (2.3) and (A.ll, since P = P, we get 



E[(u,(x2; 

2 



(2^)2 



U,(xi))(Uj(x4) - Uj(x3))] 



-2H-2„T 



ef P(k)e, {e^'^--^ 



\-2H-2„T 



e P 



(k)e, ( 



g.k-(x2-X4) _ glk-(x2-X3) _ gik-(xi-X4) 



„ik-(xi-X3; 



We use Lemma 2.2 in [4^ to get the following Fourier transform: for any x € 



(2^)2 



-2H-2i 



P(k)e*'^-"rfk = CH||x| 



\2H 



2H', 



where ch = r(l - H) / {tt2'^"+'^T{H + \)H{2H + 2)). The structure function is 
then deduced. 



Finally (2.4) coincides with the structure function obtained in [48] Section 4.5 



when u is defined by (2.2) with = oo and ^2 = 0. As this structure function 



characterizes the law of the Gaussian vector field u, this shows that the two vector 



fields defined by ( 2.3 ) and ( 2.2 ) (with ^1 = cxo and ^2 ~ 0) share the same distribution. 

□ 



A.2. Proof of Proposition 3.1. Let us denote by $0 the indicator function 



I[o,i]2(x), so that according to the construction explained in Section 13^1} the wavelets 
<f>£^s, for {l^s) G U 0, form an orthonormal basis of L^([0, 1]^). For any function 
V e L2([0, 1]2), we have for any x e [0, 1]^, v(x) = E(£,s)enuo(V' where 
(., .) denotes the scalar product in i^([0, 1]^). By the Plancherel's theorem, we deduce 
that for any k S v(k) = ^(v, $£.s)<jE>£ ^(k), where now (., .) denotes the scalar 
product in LF'i^^). Therefore, for j — 1, 2, 



(A.2) 



with 



$^,,(k)VF,(dk). 



(A.3) 



From (A.l) and the Plancherel's theorem, since the wavelets are orthogonal and nor- 
malized in L? ^ we note that 77;^ ^ are i.i.d standard Gaussian random variables. 



Now recall from (2.3) that u(x) = (ui(x), U2(x))-^ where for i = 1,2, Ui(x) 
/va(k)W^i(dk) +/v,2(k)W^2(rfk) with v,,(k) = a/(2^) (e*"^-" - l)||k||-^-i(5.j 



kikj/||k|| ). Applying (A.2) to v = Vii,Vj2, we obtain for i = 1,2 



Ui(x) 



(£,s)Gf2U0, je{l,2} 



Vij(k)$£,s(k)dk, 
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yielding the representation 



«,s)eouo 



$^.,(k)dk, (A.4) 



where 77^^^ = (vls^vlsV- 

Since the mother wavelet -0 has Af vanishing moments, for any {£, s) e ft, <f>£,s 
has Af vanishing moments along at least one direction (say xi). As a consequence, 
there exists a bounded function Oi^s such that $£^s(k) = (— ifci)*^^£^s(k) (see 
So \\m-"-^^e.si^)\ = \M~"~\-^kiV(>eA'^)\ < c||k||*^--f^^i, where c is some 
positive constant. Since A/ > iJ, the latter bound shows that k i-> ||k||^^^^$£ ^(k) 
is square-integrable on any compact. Moreover it is square-integrable at infinity as 
k I— >■ $£_s(k) is, while k 1— > ||k||^-^^^ asymptotically vanishes. Hence, for any (€, s) G 



17, k ||k|| ^ ^$£.s(k) e L^(M2) and the integral in ( |A.4p can be split, leading to 
the representation, for all x e [0, 1]^, 



u(x) = 27ra ^ [r?,,,* 



(-ff-i) 



a 

2^ 



(e*''-^ - 1) 



kk^T 



(0) 



»7ol|k 



$o(k)dk, (A.5) 



where ^\ and V are respectively defined in (3.3) and (3.4|. 



The decomposition (3.5 ) is a consequence of ( A.5 ), where e^^s — iTrar)^ g, provided 
we prove that 



(0) = 



1 



(2^)2 



R2 



(e'''-^ - 1) 



kk^ 

l|k||2 



%|lk|l-«-^$o(k)rfk. 

(A.6) 



Since by assumption Jj^ j^j2 u(x)(ix = 0, integrating both sides of (A.5 1 on [0, 1]^ 
leads to 



E / ^ 

1 



i-H-l) 



(x)dx - V 



i-H-l) 



(0) 



dx / (e*''-^ - 1 



kk 



|k||2 



77ol|k|r«-^$o(k)dk = 0. (A.7) 



Let u;(x) = 17^ g^^^ (^)- Since wavelets possess at least one vanishing moment 



/jg ^2 If (x)dx = that is to say w{Q) ~ 0. According to definition (3.4) of the Leray 
projector this implies that Jr„ Vw{n)dx. — and therefore 



E 

(£,s)6n 



[0,1]^ 



i-H-l) 



(x)dx = 0. 



(A.; 
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Now consider the last term in (A. 7). We have for j — 1,2 
(e*"^-^ - 1) 



d 



r,o||k|r^-i<i>o(k)dk 



kk^ 



r7o||k|r^-i$o(k)dk 



kk^ 



rjollkll 



-H-1 



-ikx 



d_ 
dxi 



I[o^i]2 (x)(ix j dk = 



Therefore, we obtain 



[0,1]2 



dx / (e^'^ '' - 1) 



kk^ 



r7o|lk|r«-^$o(k)dk 



(e^*"" - 1) 



kk 



T 1 



IkiP 



J7ol|k|r^-i$o(k)dk. 



Using (A. 8) and (A.9) in (A.7) proves (A.6), which concludes the proof. 



(A.9) 



□ 



A. 3. Proof of proposition 3.2, From (3.5) and the definition of V (see (3.4)), 
we have 



Since *^,s e iL([0. ^e have V [*, 



*£.s- So 



Since e(ij) are iid zero-mean Gaussian random variables with variance (27rCT)^I, we 
have 



Let us simplify the sum above. First, recall that $^ ^ "^(x) = (—A) 2 $^ s(x), so 
for any fc = 1, 2: 



= ( E *ij(*iJ,(-A)"^^^\,,),(-A)^§,%). 
Since the mother wavelet iJj ha s Af > vanishing moments, similar arguments as 



in the proof of proposition 



3.1 



lead to |||k||^^-i*^-, ^,(k)| < c||k 



\M-H-1 



where c is 
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some positive constant. So /jpy2(-A) 2 v[/fe ^, (x)dx = ||k|| ^vj/^, g,(k) =0 
and 

(ij)60 (i,j)Gnuo 
Therefore 

(ij)ef2 

and 

E[d,,,d,.,,,] = (2^a)2((-A)^*,,„ (-A)"^ (A.IO) 
Second, let us introduce the curl operator, defined for any / e L^{[0, 1]^) by 

curl f{xi,X2) ^ i- J f{xi,t)dt, J f{t,X2)dt 

and note that ^£.s = curl[<I>^_s]- It is easily checked that the adjoint operator is 
given for g = (.91,52)^ by curl g{xi,X2) = --£^gi{xi,X2) + gf^ 52(2^1, 2^2) and that 
curl curl = (— A)^^. Consequently we have from (A.IO) 

= (2^a)2(cTirl*c'iirl[$^7/-')], $^rif"'^) 

= (2..)^((-A)-ci>(/-),<i>(7r^)) 

= (2..)^((-A)-i$(/-),(-A)-^$(7--)) 

where the existence of wavelet $^ ^ is guaranteed by the M > H vanishing 
moments of the mother wavelet ip. 

□ 



A. 4. Proof of proposition 3.3 , By definition of S and S ^, and by use of the 

biorthogonality of the wavelets family {^^/^s^'^^', (•^' ^) ^ ^1' have for any a G £^($7) 

(f',s')eo (ij)eo 
(^',s')eo 

We can show similarly that S^^Sa = a. Therefore operator corresponds to the 
inverse operator of S. 

□ 
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A. 5. Proof of proposition 5.1 The gradient with respect to ^ of the data- 
term (5?/(e„) in (4.121 is given by inner products with the fractional divergence-free 



wavelets. Indeed, we have: 



de] 5y{en) = ((yi(x, e„) - 2/o(x)), Vj/i(x, e^f 



d 



d 

(yi(x, e„) -yo(x))Vyi(x, e„)/7rT— ■P 



■V 



i-H-l) 



(x)) 
(x)> 



and by Fourier-Plancherel formula: 



d 



dej jyi^n) = ((yi(k,e„) - j/o(k))V2/i(k,e„)/^ 



e^.sllk 



= ((yi(k,e„)-2/o(k))Vyi(k,e„)/l|k 



||fe||2 
fclfc2 

ll*:|P , 



Q,s$£,s(k)) 



fclfc2 



(yi(k, e„) - yo(k))Vyi(k, e„), $^,s(k)), 



= ((-A)"^7'i[(yi(-, e„) - yo(-))Vyi(-, e„)](x), $,,,(x)). 

The gradient with respect to ^ is obtain similarly. The gradient of the regularization 
term is simply: 9e„7^(e„) = ji^^^^n- □ 



A. 6. Proof of proposition 5.2, The gradient (5.51 of the data-term Sy{dn) is 
obtained analogously to (5.1 1. For the gradient of the regularizer term, by definition 
Q of S-^: 



1 



/3(27r(T)2 



Formula (5.6 1 follows from definition of see (5.4). 



□ 



Appendix B. Adaptation of algorithm [3] for irregular wavelets. 

If tp is not 2H + 4 but only H + 2 times differentiable, one may replace iii)-v) in 
algorithm [3] by the following steps: 

- compute A^($^'''*'-'(i„) in the Fourier domain by FFT 

- compute the FWT using orthogonal wavelets s) G 17„} to get the 
nxn matrix denoted by |e„] whose element at row index si) and column 
index (£2,^2) is (A^(*i'(°)d„), 

- compute off-line matrices F^") for a = 0, 1,2 with algorithm [H] (section 5.2.21 

- obtain the scalar product in (|5.6|) by addition of matrix products 



l,(H+2) 



(e'.s')en„ 



= (F(2)"[e„lF(°) -f 2F(i)"[e„lF(i) +F(0)"le„lF(2)) 
where (M)^.^ denotes the (£, s)-th element of matrix M 
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Appendix C. Matrices of mono-dimensional connection coefficients. 

The matrices F*^") involved in (5.121, where < a < + 2, are composed of 



wavelets connection coefficients defined in (5.8 1. Note that F'"' is the identity matrix 



since we are considering an orthonormal basis. Moreover, for a being a positive 
integer, fractional Laplacian operator becomes a standard differentiation up to factor 
(— 1)" and F^") can be computed by solving an eigenvalue problem as detailed in 
[5J[T2]. However, in the more general case of fractional Laplacian differentiation, no 
method have been explicitly proposed in literature. In this appendix we provide an 
approximation of F^"^ in terms of scaling functions connection coefficients, that turn 
out to be easily computable as the solution of a linear system, as explained in the 
following. 

C.l. Matrix F^"^. In this section we assume a > and we show that any 
entry /|"''^, of F*^") can be determined recursively from an infinite series of connec- 
tion coefficients of scaling functions defined at the finest scale s„ = log2{n). These 
connection coefiicients, denoted by ei"', are given for any i?, e Z by 



(C.l) 



where (f denotes the scaling function associated to wavel et ^p. An efficient algorithm 
for the computation of ei"' is obtained in section C.2 We hereafter propose an 



approximation of f^"\, g, as a truncation of these infinite series of scaling function 
connection coefficients. 

Let us begin by recalling the two-scale relations associated to the orthonormal 



wavelet basis {tpe,s{x);x e E,£, s S Z} defined in (3.1), see 

(p{x) = \/2^ hk(p{2x - fc), 

k 

ijjix) = V2^^gkip{2x - k), 



(C.2) 
(C.3) 



where h and g are the conjugate mirror filters of finite impulse response given by hk = 
^{ip{x),Lp{2x-k)) andfffc = -^{ip{x),Lp{2x-k)). For any function &(£, f) e L'^{I?), 



let us define the following convolution operators: 

Gi&(^,/) ^^.9fc6(2^ + fc,f). 

k 

GX£,£')=Y.9kb{e,2£' + k) 

k 

H^h{£,£') = V2^hkb{2£ + k,£') 

k 

H^h{£,£') ^ V2^hkb{£,2£' + k). 



(C.4) 
(C.5) 
(C.6) 
(C.7) 



We also consider operator h}^^ (resp. I]j2 ^) defined by iterating i times operator iJj^ 
(resp. H_2)- Following the methodology introduced in [7] (see details in [39]), we 
obtain from (C.2)-(C.7) that for s, s' < s„ 



2 \ " 



(C.8) 
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To get a similar representation for fi"\, ^, , we need to consider the same procedure 
with periodized wavelets and scaling functions instead of ^ and ip. It can be shown 
that in the case of scaling functions defined at scale s„ and periodized over [0,1], 
connection coefficients are (2'*")-periodic functions 



+00 

E 

fc.A:' — — 00 



e + k), 



+00 



fc— — 00 



provided the latter series converges. Therefore, by redefining operators (C.4)-(C.7| 
with circular convolution on (2*")-periodic signals, we obtain similarly as (C.8): for 
< s,s' < s„ and for < ^<2''-i,0 < l'<2^'-'^, 



Aa) — r n 



±J-2 



+00 
A: — —00 



(C.9) 



provided the latter series is convergent. The remaining terms of F^"^ can be treated 
in the same way: for < s < s„, < £<2'*~^ 



Jo,o,e,s 



Ac,) 



Ac) 

^0,0,0,0 



_2i3-l il.2 



-\-oo 

E ' 

fc— — 00 

+00 

A; — — 00 



+00 



k — — CK 



s(")(fc2^",0) 



(C.IO) 



Recursive formulae (|C.9|)-(|C.10P show that the knowledge of e 
the matrix F 

" Ac) 



(a) 



entirely determines 



Finally, as it will be explained in section 



where Ca > 0. Since ei"^(£,^') — e 

andforany fc 9^ 0, ei"^(£+A:2"",f ) behaves as c„|£'-^-/c2" 
large. This shows that the terms associated to fc ^ in ( C.9 )-( C.IO ) are negligible with 
respect to the the terms associated to fc = 0, provided 2^" is sufficiently large. The 
latter condition is a reasonable assumption in standard image setting where typically 
Sra > 8. This is the reason why we propose the following approximation, for any 
< £,f <2''": 



- J") I 



C.2 



(^,o)^c„K|-i- 



as £ 



f ,0), we deduce that for any < <2^", 
-2a j£ jg sufficiently 



+00 

E 

k— — oo 



• + k2'",e') 



Ac) 



for 0<\e -£\< 2" 



es„ (^:^' - 2"") for 2''"-^ < \i' - £\ < 2^" 



(C.ll) 



This approximation is based on the above explanation when \£' — i\ < 2" 
extended to 2'^"~^ < \i' — £\ < 2*" in order to respect (2''")-periodicity. 



and is 



The derivation of matrix F*^"^ is thus very simple: from (C.9|-(C.10), we see 
that matrix F^"-* is a bi-dimensional anisotropic discrete wavelet transform of the 



(2^")-periodic function X^fc^oo + ^2*",^'), where the latter is approximated 

by (C.ll). In other words, relations (C.9|-(C.10) perform a basis change from the 
orthonormal family {Efe fe' ¥'(2"" (xi + k) - e)ip{2''^{x2 + k) - £');0 < £,£' < 2""} 
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to the orthonormal family {<i>^^s; (£, s) G iln U 0}. Indeed, recursive convolutions 
appearing in (C.9|-(C.10) implement (up to the multiplicative constant 2^*") the fast 
recursive filtering algorithm proposed by Mallat [531 for FWT. In practice, we thus 
compute F^"' by a simple FWT of the discrete function defined in the right hand side 
of dUHj ). 

C.2. Computation of connection coefficients ei"\ We hereafter adapt the 
general framework proposed by Beylkin in [5, 6J to the case of the computation of 
scaling function connection coefficients appearing in (C.II). 

The fractional Laplacian operator is rewritten as a convolution operator for any 
scaling function with a compact support. Indeed, if a — f /2 € M\N, fractional Lapla- 
cian can also be defined by Riesz potentiaQ [20J : 



with Cn 



v/^r(-a)2-^° 
r((l+2a)/2) 



ip(x) 



+ 00 



(p{z) 



\x-z\^ 



-dz. 



In the previous expression, the convolution kernel writes 



k{x) 



1 



Ca\x\ 



2a+l ' 



(C.12) 



Since we have Cs"' = ei"\£ — £',0), the computation of all scaling function 
connection coefficients reduces to the computation of ei"\£, 0) for £ gZ. From (C.2), 
we derive that 



^^^\x) = V2j2hk 



fc=0 



-92 
dx"^ 



L-l 



Lp{2x - fc) = V22^°' hkip^'^\2x - k), 



k=Q 



where L denotes the number of non zero coefficients of the scaling filter h. Using 
(C.2 1 for if and the above relation for tp^°''> in (C.l) leads to 



L-l L-l+2e-k 



ei:)(£,0) = 22"^ ^ hkhk-2i+jeizHj,0). 

k=0 j=2i~k 



k-2i+je):;'(j,u). (C.I3) 
Moreover, an asymptotic behavior can be derived from the Taylor expansion of 

1 



the kernel ( [CT2l ) as in [6]: for £ ^ cx), 

1 



\£\l+2a~ 



-2M 



(C.14) 



In order to compute ei"''(£, 0), we solve the linear system (C.13) subjected to the 
above asymptotic behavior as boundary conditions. Specifically, for \£\ > £mim where 
£min is chosen sufficiently large (typically ^min > n/8), we set ei"\£,0) = ^ |£|\+2q ■ 

Then for £ — — Anin, •■•,^min, an analytical solution ei'^\£,0) of (C.13) is obtained as 
described below. 

Let Hk be the function defined for any of £,j £ Z by 



Hk{£,j) 



hkhk-2e+3 if 2£-k<j<L-l + 2£-k, 
otherwise. 



*This definition can be extend toa — 1/2GN using some appropriate kernel, see e.g. |42| 
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Let Hfe be the matrix of size {2£min + 1) x i^^min + 1) whose element at row I and 
column j is H^it — imin^j — ^min)- Let I denote the identity matrix. Then 



e(^) = [22" ^ Hfc - I]-ib (C.15) 



fe=0 



where e^"^ and b are (2£niin+l)-dimensional vector whose components are respectively 

ei:'(£,0) and 

~ ~ 2^ 2^ „ |„-|l+2a 2^ 2^ p h-|l+2a 
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Fig. C.l. Synthesized fBms u(x) = (tii(x), U2(x^ 
vorticity maps, i.e. 9a;ji2(x) — dyUi{x) (right) . 
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Fig. C.3. Performance of the regularizers according to the value of H . Regularization coefficient 
for methods A to E (see section \6.S\ for a description) were chosen to minimize the RMSE. The 
given criteria are in order RMSE/MBAE/SAE. For each value of H , the 3 best results with respect 
to each criterion are displayed in blue. 
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Fig. C.4. Radial power spectrum estimates in logarit hmic coordinates (above) and ordinary 
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least squares fitting (below) for methods A to E (see section 
and H = 1 (right), compared to ground truth and the theoretical decay, 
where chosen to minimize the RMSE. 
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Fig. C.5. Estimated vorticity maps for methods A to E (see section \6.3\ for a description) and 
for J/ = I (above) and H = 1 (below). Regularization coefficients were chosen to minimize the 
RMSE 
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Fig. C.6. Influence of regularization coefficient on the accuracy of the estimate in terms of 
RMSE (above), MBAE (middle) and SAE (below) for methods A to E (see section [gT^] for a de- 
scription) and for H = ^ (l^ft) cifid, H = 1 (right). 



